/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
/** @file Lanczos.h
 *  @brief Definition of class Lanczos - solution of the eigenvalue problem
 *
 *  EIGENVALUE PROBLEMS:  Class LANCZOS  
 *  STANDARD AND GENERALIZED EIGENVALUE PROBLEMS 
 *  A*x = L*x and A*x = L*B*x
 * 
 *  Lanczos transformation 
 */ 

#include "Lanczos.h"
namespace smv_lib
{

//! empty (default) constructor
Lanczos::Lanczos(void)
{
  isMatrix_A   = false;
  isMatrix_B   = false;
  isCalculated = false;
}

//! constructor for the standard problem
Lanczos::Lanczos(Matrix& a, int n) : A(a), No_Lanc(n)
{
  N = A.getRows();
  assert(N == A.getCols());
  assert(No_Lanc <= N);
  isMatrix_A   = true;
  isMatrix_B   = false;
  isCalculated = false;
}

//! constructor for the generalized problem
Lanczos::Lanczos(Matrix& a, Matrix& b, int n) : A(a), B(b), No_Lanc(n)
{
  N = A.getRows();
  assert(N == A.getCols());
  assert(No_Lanc <= N);
  int N_B = B.getRows();
  assert(n == B.getCols() && N == N_B);
  isMatrix_A   = true;
  isMatrix_B   = true;
  isCalculated = false;
}

//! destructor
Lanczos::~Lanczos(void)
{}

//! set Matrix A (for the standard problem)
void Lanczos::setMatrix_A(Matrix& a)
{
  A = a;
  N = a.getRows();
  assert(N == a.getCols());
  isMatrix_A   = true;
  isMatrix_B   = true;
  isCalculated = false;
}

//! set Matrices A and B (for the generalized problem)
void Lanczos::setMatrix_B(Matrix& a, Matrix& b)
{
  A = a;
  B = b;
  N = A.getRows();
  assert(N == A.getCols());
  int m = B.getRows();
  assert(m == B.getCols() && N == m);
  isMatrix_A   = true;
  isMatrix_B   = true;
  isCalculated = false;
}

//! set number of eigenvalue/vector pairs
void Lanczos::setNumber(int n)
{
  No_Lanc = n;
  assert(No_Lanc <= N);
}

//! check and calculate
void Lanczos::calculate()
{
  Matrix Modal(N,No_Lanc,"Lanc_Modal.txt");
  Vector Spektar(No_Lanc,"Lanc_Spec.txt");
  V = Modal;
  d = Spektar;
  //--- standard or generalized problem
  if (isMatrix_A == true && isMatrix_B == false)
    stLanczos();
  else if (isMatrix_A == true && isMatrix_B == true)
    genLanczos();
  //--- sort eigenvectors and eigenvalues
  eigsrt_Lanc();
  isCalculated = true;

  //--- write obtained results to files
  V.Write();
  d.Write();
  cout << "\nEigen-Vectros are in the file LANC_MODAL.txt\n";
  cout << "\nEigen-Values are in the file LANC_SPEC.txt\n";
}

//! give the modal matrix
Matrix& Lanczos::giveModalMatrix()
{
  Matrix* Zero = new Matrix(N,N);
  if(isCalculated == true)
    return V;
  else
    return (*Zero);
}

//! give eigenvalues
Vector& Lanczos::giveSpectrum()
{
  Vector* zero = new Vector(N);
  if(isCalculated == true)
    return d;
  else
    return (*zero);
}

//-------- Standard Lanczos Transformation
//------------------------------------------------------------------------------
/**
** \brief stLanczos - Lanczos transformation for the standard eigenvalue problem
**
**  Transformation of Matrix A(n,n) into tridiagonal form, 
**  using the Lanczos algorithm
**
**/
//------------------------------------------------------------------------------
void Lanczos::stLanczos()
{
  int m;
  m = (int)ceil(1.5 * No_Lanc);
  cout << "\n----- m = " << m << endl;
  standLanczos(m); // standard Lanczos with No_Lanc vectors
}

//---- standard Lanczos with m Lanczos vectors
//
void Lanczos::standLanczos(int m)
{
  int i,j;
  Matrix* Tm = new Matrix(m,m);
  assert(Tm != NULL);
  Tm->Identity();
  Vector rj(N);                 // intermediate vector rj (with N elements equal to zero)
  //---- vectors r
  Vector* r = new Vector[m+1];  // array of temporary vectors
  assert(r != NULL);
  //--- initialize vectors
  for (i=0; i<m+1; i++) {
    r[i].setSize(N);
    for (j=0; j<N; j++) {
      if (i==0)
        r[i].setData(j, 1.0);  // initial vector r_0 with all elements equal to 1.0
      else
        r[i].setData(j, 0.0);
    }
  }
  //---- Lanczos vectors q
  Vector* q = new Vector[m+1];  // array of m Lanczos vectors
  assert(q != NULL);
  for (i=0; i<m+1; i++) {  // all Lanczos vector initialized to zero
    q[i].setSize(N);
    for (j=0; j<N; j++) {
      q[i].setData(j, 0.0);
    }
  }
  double* Alpha = new double[m+1];   // diagonal elements of matrix Tm
  assert(Alpha != NULL);
  double* Beta  = new double[m+1];  // off-diagonal elements of matrix Tm
  assert(Beta != NULL);

  //--- start the process
  double factor;
  //cout << "\nMatrix A = \n" << A << endl;
  Alpha[0] = 0.0;
  Beta[0] = 0.0;
  //============================================== initial orthonormalization
  Beta[1] = r[0].Norm();
  factor = 1.0 / Beta[1];
  Vector temp1(N), temp2(N), temp3(N);
  temp1 = r[0].Dot(factor);
  //=============================================== the first Lanczos vector
  for (i=0; i<N; i++) {
    q[1].setData(i,temp1.getData(i));
  }

  //=============================================== other Lanczos vectors
  for (j=1; j<m; j++) {   // (j=1,2,...,m-1)
    //=================================================== solution of eqs.
    //rj = A_inv * q[j];
    cout << "j = " << j << "--- rj = " << rj << endl;
    //=================================================== orthogonalization
    //Alpha[j] = q[j] * rj;
    Alpha[j] = q[j].Dot(rj);
    cout << "\n--- Alpha[j] = " << Alpha[j] << endl;
    //Beta[j]  = r[j-1].Norm();
    cout << "\n--- Beta[j]  = " << Beta[j] << endl;       // ok
    temp1 = q[j].Dot(Alpha[j]);    // q_j * alfa_j
    temp2 = q[j-1].Dot(Beta[j]);   // q_(j-1) * beta_j
    temp3 = temp1 + temp2;
    for (i=0; i<N; i++)
      r[j].setData(i, r[j].getData(i) - temp3.getData(i));
    //r[j] = rj - temp3;
    //r[j] = rj - q[j].Dot(Alpha[j]) - q[j-1].Dot(Beta[j]);
    cout << "\n--- r[j]     = " << r[j] << endl;        
    //=================================================== orthonorming i Lanczos vectors
    Beta[j+1] = r[j].Norm();
    factor = 1.0 / Beta[j+1];
    temp1 = r[j].Dot(factor);
    for (i=0; i<N; i++)
      q[j+1].setData(i,temp1.getData(i));
    //q[j+1] = r[j].Dot(factor);
  }
  //=================================================== the last Alpha coeff.
  //rj = A_inv * q[m];
  //Alpha[m] = q[m] * rj;
  Alpha[m] = q[m].Dot(rj);

  cout << "\n----------Standard Lanczos\n" << endl;
  for (i=0; i<m+1; i++)
    cout << "-- Alpha[" << i << "] = " << Alpha[i] << endl;
  for (i=0; i<m+1; i++)
    cout << "--  Beta[" << i << "] = " << Beta[i] << endl;
  cout << "\n--------- Lanczos vectors\n" << endl;
  for (i=0; i<m+1; i++)
    cout << " i = " << i << "  => q[" << i <<"] = " << q[i] << endl;

}

void Lanczos::genLanczos()
{}
//================ eigsrt
//------------------------------------------------------------------------------
/**
**  \brief eigsrt: sorting eigenvalues/eigenvectors in ascending order
**
**  INPUT:
**  Unsorted eigenvalues in Vector Lam(n) and
**  corresponding eigenvector matrix X(n,n)
**  (as obtained from JACOBI or TQLI)
**  OUTPUT:
**  Sorted into ascending order eigenvalues in Vector Lam(n) and
**  correspondingly sorted eigenvector matrix X(n,n)
**  Note: Vector Lam(n) ans Matrix X(n,n) are the private members of EigVal class
**/
//-------------------------------------------------------------------------------
void Lanczos::eigsrt_Lanc()
{
  int k, j, i;
  double p;
  int n = d.getLen();

  for (i=0; i<n-1; i++)  {
    k = i;
    p = d.getData(k);               // p = d(k)
    for (j=i+1; j<n; j++) {
      if (d.getData(j) <= p) {    // d(j) <= p
        k = j;
        p = d.getData(k);       // p = d(k)
      }
    }   // end for j = i+1 ...
    if (k != i)  {
      d.setData(k,d.getData(i));   // d(k) = d(i)
      d.setData(i,p);                // d(i) = p
      for (j=0; j<n; j++)  {
        p = V.getElement(j,i);                  // p = V(j,i)
        V.setElement(j,i,V.getElement(j,k));    // V(j,i) = V(j,k)
        V.setElement(j,k,p);                    // V(j,k) = p
      }
    }  // end for k != i
  }   // end for i = 0 ...
}


}   // end of namespace smv_lib
//-------------------------------------- end of file Lanczos.cpp
